#安装
{ install.packages("purrr")
  install.packages("furrr")
  install.packages("devtools")
  BiocManager::install("ComplexHeatmap")
  BiocManager::install("igraph")
  if(!require(dplyr))install.packages("dplyr")
  BiocManager::install(c("SingleR", "celldex"))
  BiocManager::install("scRNAseq")
  devtools::install_github("jinworks/CellChat")
  BiocManager::install(c(
    "pathview",       # 用于KEGG通路可视化
    "ComplexHeatmap", # 用于高级热图绘制
    "circlize",       # 提供colorRamp2函数
    "ggraph",         # 用于网络可视化
    "igraph"          # 网络分析基础包
  ))
  install.packages(c(
    "ggplot2",
    "dplyr",
    "tidyr"
  ))
  # 安装Harmony（如果尚未安装）
  if (!require("harmony"))  install.packages("harmony")
}
#加载
{
  library(harmony)
  library(pathview)
  library(circlize)
  library(ggraph)
  library(tidyr)
  library(Seurat)
  library(SingleR)
  library(celldex)
  library(scRNAseq)
  library(dplyr)
  library(ggplot2)
  library(EnhancedVolcano)
  library(ggrepel)
  library(ComplexHeatmap)
  library(igraph)
  library(CellChat)
  library(patchwork)
  library(harmony)
  library(purrr)
}
setwd("C:/Users/13359/Desktop/code&output/数据结果/DEG/")
# 基本分析
{
  # 设置主路径（根据你的实际路径修改）
  main_dir <- "E:/R-date/GSE190510_RAW"
  
  # 1. 遍历子文件夹并创建带过滤的Seurat对象
  subdirs <- list.dirs(main_dir, full.names = TRUE, recursive = FALSE)
  seurat_list <- list()
  
  for (dir_i in subdirs) {
    # 读取每个样本的10X数据
    mat <- Read10X(data.dir = dir_i)
    
    # 创建带初步过滤的Seurat对象
    sample_name <- basename(dir_i)
    seurat_list[[sample_name]] <- CreateSeuratObject(
      counts = mat,
      project = sample_name,
      min.cells = 3,      # 基因至少在3个细胞中出现
      min.features = 200  # 细胞至少检测到200个基因
    )
    
    # 添加线粒体基因比例
    seurat_list[[sample_name]][["percent.mt"]] <- PercentageFeatureSet(
      seurat_list[[sample_name]], 
      pattern = "^MT-"
    )
  }
  
  # 2. 合并所有样本
  combined_seurat <- merge(seurat_list[[1]], y = seurat_list[-1])
  
  # 3. 执行细胞过滤
  combined_seurat <- subset(
    combined_seurat,
    subset = nFeature_RNA > 500 &
      nFeature_RNA < 3000 &
      percent.mt < 10 &
      nCount_RNA > 500
  )
  combined_seurat <- JoinLayers(combined_seurat, layers = "counts")
  counts_matrix <- GetAssayData(combined_seurat, slot= "counts")
  gene_counts <- Matrix::rowSums(counts_matrix > 0)
  total_cells <- ncol(counts_matrix)
  threshold <- 0.001 * total_cells  # 计算0.1%阈值
  genes_to_keep <- names(which(gene_counts >= threshold))
  combined_seurat <- subset(combined_seurat, features = genes_to_keep)
  
  # 查看最终元数据
  head(combined_seurat@meta.data, 12)
  VlnPlot(combined_seurat,features = c("nFeature_RNA","nCount_RNA","percent.mt"),ncol = 3)
  plot1 <- FeatureScatter(combined_seurat, feature1 = "nCount_RNA", feature2 = "percent.mt")
  plot2 <- FeatureScatter(combined_seurat, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
  CombinePlots(plots = list(plot1, plot2))
  combined_seurat <- JoinLayers(combined_seurat, layers = "counts")
  # 提取全局counts矩阵
  counts_matrix <- GetAssayData(combined_seurat, layer = "counts")
  combined_seurat <- NormalizeData(combined_seurat,normalization.method = "LogNormalize",scale.factor = 10000)
  combined_seurat <-FindVariableFeatures(combined_seurat,selection.method = "vst",nfeatures = 3000)
  top10 <-head(VariableFeatures(combined_seurat),10)
  plot1 <-VariableFeaturePlot(combined_seurat)
  plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)  
  plot1+plot2
  combined_seurat <- FindVariableFeatures(combined_seurat, selection.method = "vst",nfeatures = 3000)
  combined_seurat <- ScaleData(combined_seurat)
  combined_seurat <- ScaleData(combined_seurat, features = VariableFeatures(combined_seurat))
  combined_seurat <-RunPCA(combined_seurat,features = VariableFeatures(object = combined_seurat))  
  print(combined_seurat[["pca"]],dims=1:5,nfeatures=5)
  VizDimLoadings(combined_seurat,dims = 1:2,reduction = "pca")
  DimHeatmap(combined_seurat,dims = 1,cells = 500,balanced = TRUE)
  DimHeatmap(combined_seurat,dims = 1:10,cells=500,balanced = TRUE)
  combined_seurat<- JackStraw(combined_seurat, dims = 20, num.replicate = 100)
  combined_seurat <-ScoreJackStraw(combined_seurat,dims = 1:20)
  JackStrawPlot(combined_seurat,dims = 1:20)
  ElbowPlot(combined_seurat)
  #harmony前的去批次所得步骤
  combined_seurat <-FindNeighbors(combined_seurat,dims = 1:10)
  combined_seurat <-FindClusters(combined_seurat,resolution = 0.4)
  combined_seurat <-RunUMAP(combined_seurat,dims = 1:10)
  DimPlot(combined_seurat,reduction = "umap")
  #去批次后umap
  print(table(combined_seurat@meta.data$orig.ident))
  combined_seurat <- RunHarmony(
    combined_seurat,
    group.by.vars = "orig.ident", 
    assay.use = "RNA",
    theta = 2,
    lambda = 1,
    max.iter.harmony = 20
  )
  # 2. 使用Harmony降维结果重新进行聚类和可视化
  combined_seurat <- FindNeighbors(combined_seurat,reduction = "harmony", dims = 1:10)
  combined_seurat <- FindClusters(combined_seurat,resolution = 0.4)
  combined_seurat <- RunUMAP(combined_seurat,reduction = "harmony",dims = 1:10,reduction.name = "umap_harmony")
  DimPlot(combined_seurat,reduction = "umap_harmony")
  combined_seurat <-RunTSNE(combined_seurat,reduction = "harmony",dims=1:10,reduction.name = "tsne_harmony")
  DimPlot(combined_seurat,reduction =  "tsne_harmony")
  Layers(combined_seurat[["RNA"]])
  combined_seurat <- JoinLayers(object = combined_seurat, layers = "data")
}
save.image("基本分析.RData")
saveRDS(combined_seurat,file = "基本分析seurat.rds")
#细胞注释
{
    ref_data <- celldex::HumanPrimaryCellAtlasData()  # 通用人类细胞类型
    # 小鼠参考数据集（若为小鼠肝脏）
    # ref_data <- celldex::MouseRNAseqData()
    expression_matrix <- GetAssayData(combined_seurat, layer = "data")  # 必须用归一化后的数据
    # 提取参考集的标签
    labels <- ref_data$label.main  # 粗粒度注释
    # 运行SingleR
    singleR_results <- SingleR(
      test = expression_matrix, 
      ref = ref_data, 
      labels = labels,
      assay.type.test = 1  # 确保使用log归一化数据
    )
    # 添加注释到metadata
    combined_seurat$SingleR_labels <- singleR_results$labels
    
    # 检查前6个细胞的注释结果
    head(combined_seurat@meta.data$SingleR_labels)
    table(combined_seurat@meta.data$SingleR_labels)
    library(ggplot2)
    p_umap <- DimPlot(combined_seurat, 
                      reduction = "umap_harmony",
                      group.by = "SingleR_labels",
                      label = TRUE, 
                      pt.size = 0.5) + 
      ggtitle("SingleR Cell Annotation (UMAP)")
    ggsave("SingleR_UMAP.pdf", p_umap, width = 12, height = 8)
    p_tsne <- DimPlot(combined_seurat, 
                      reduction = "tsne_harmony",
                      group.by = "SingleR_labels",
                      label = TRUE,
                      pt.size = 0.5) + 
      ggtitle("SingleR Cell Annotation (t-SNE)")
    ggsave("SingleR_tSNE.pdf", p_tsne, width = 12, height = 8)
  }
save.image("细胞注释.RData")
saveRDS(combined_seurat,file = "细胞注释seurat.rds")
load("细胞注释.RData")
# 添加分组标签
{
  unique_ids <- unique(combined_seurat$orig.ident)
  combined_seurat$group <- ifelse(
    combined_seurat$orig.ident %in% unique_ids[1:5], 
    "ASS-ILD", 
    "Healthy"
  )
}
#细胞类型的分布
{
  # 统计各组细胞类型数量
  celltype_counts <- table(combined_seurat$SingleR_labels, combined_seurat$group)
  print("原始计数表:")
  print(celltype_counts)
  
  # 计算比例（按列计算，即每组内部的比例）
  celltype_props <- prop.table(celltype_counts, margin = 2) * 100
  print("比例表(%):")
  print(round(celltype_props, 2))
  
  # 转换为数据框便于ggplot绘图
  celltype_df <- as.data.frame(celltype_counts)
  colnames(celltype_df) <- c("CellType", "Group", "Count")
  
  # 添加比例列
  prop_df <- as.data.frame(celltype_props)
  colnames(prop_df) <- c("CellType", "Group", "Percentage")
  celltype_df <- merge(celltype_df, prop_df, by = c("CellType", "Group"))
  library(ggplot2)
  library(RColorBrewer)
  
  # 设置颜色
  n_celltypes <- length(unique(celltype_df$CellType))
  celltype_colors <- colorRampPalette(brewer.pal(8, "Set2"))(n_celltypes)
  
  # 堆叠条形图（计数）
  p_count <- ggplot(celltype_df, aes(x = Group, y = Count, fill = CellType)) +
    geom_bar(stat = "identity", position = "stack") +
    scale_fill_manual(values = celltype_colors) +
    labs(title = "Cell Type Distribution by Group (Count)",
         x = "Group", y = "Number of Cells") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1),
          plot.title = element_text(hjust = 0.5))
  
  # 堆叠条形图（百分比）
  p_percent <- ggplot(celltype_df, aes(x = Group, y = Percentage, fill = CellType)) +
    geom_bar(stat = "identity", position = "stack") +
    scale_fill_manual(values = celltype_colors) +
    labs(title = "Cell Type Distribution by Group (Percentage)",
         x = "Group", y = "Percentage (%)") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1),
          plot.title = element_text(hjust = 0.5))
  
  # 并排显示
  library(patchwork)
  combined_plot <- (p_count | p_percent) + plot_layout(guides = "collect")
  print(combined_plot)
  ggsave("CellType_Distribution_by_Group.pdf", combined_plot, width = 14, height = 6)
}
#组间差异分析
{
  # 创建存储目录
  if(!dir.exists("DEG_results")) dir.create("DEG_results")
  
  # 设置待分析的细胞类型 (注意名称需与metadata中完全一致)
  target_celltypes <- c("T_cells", "NK_cell", "B_cell", "Monocyte")
  
  # 创建组合标识并设置细胞身份
  combined_seurat$celltype.group <- paste(
    combined_seurat$SingleR_labels,
    combined_seurat$group, 
    sep = "_"
  )
  Idents(combined_seurat) <- "celltype.group"
  
  # 差异分析函数 (保持不变)
  run_deg_analysis <- function(seurat_obj, cell_type, group1 = "ASS-ILD", group2 = "Healthy") {
    ident1 <- paste(cell_type, group1, sep = "_")
    ident2 <- paste(cell_type, group2, sep = "_")
    
    deg_results <- FindMarkers(
      object = seurat_obj,
      ident.1 = ident1,
      ident.2 = ident2,
      test.use = "wilcox",
      min.pct = 0.25,
      logfc.threshold = 0.25,
      only.pos = FALSE,
      verbose = FALSE
    )
    
    deg_results$gene <- rownames(deg_results)
    deg_results$comparison <- paste(group1, "vs", group2)
    deg_results$cell_type <- cell_type
    return(deg_results)
  }
  
  # 只分析目标细胞类型
  all_deg_results <- list()
  for(cell in target_celltypes) {  # 修改为遍历指定细胞类型
    deg_df <- tryCatch(
      run_deg_analysis(combined_seurat, cell),
      error = function(e) {
        message(paste("Error in", cell, ":", e$message))
        NULL
      }
    )
    
    if(!is.null(deg_df)) {
      write.csv(deg_df, paste0("DEG_results/", cell, "_ASS-ILD_vs_Healthy.csv"), row.names = FALSE)
      all_deg_results[[cell]] <- deg_df  # 使用细胞类型作为列表名
    }
  }
  
  # 合并结果时处理空值
  combined_deg <- do.call(rbind, compact(all_deg_results))  # 使用purrr::compact过滤空值
  
  
  # 修改后的可视化函数
  plot_volcano <- function(deg_df, title) {
    # 计算上下调基因数量（p_adj < 0.01且|logFC| > 0.5）
    up_genes <- sum(deg_df$p_val_adj < 0.01 & deg_df$avg_log2FC > 0.5, na.rm = TRUE)
    down_genes <- sum(deg_df$p_val_adj < 0.01 & deg_df$avg_log2FC < -0.5, na.rm = TRUE)
    
    # 在标题中添加统计信息
    new_title <- paste0(title, "\nUp: ", up_genes, " genes | Down: ", down_genes, " genes")
    
    EnhancedVolcano(
      deg_df,
      lab = deg_df$gene,
      x = "avg_log2FC",
      y = "p_val_adj",
      pCutoff = 0.05,
      FCcutoff = 0.5,
      title = new_title,
      subtitle = "",
      caption = NULL
    )
  }
  
  # 修改后的点图函数
    
    # 修改后的点图函数 - 仅展示差异最大的前30个基因
  plot_dotplot_top_genes <- function(seurat_obj, cell_type, n_genes = 30) {
    subset_data <- subset(seurat_obj, SingleR_labels == cell_type)
    
    cell_deg <- combined_deg %>% 
      filter(
        cell_type == !!cell_type,
        p_val_adj < 0.01,          # 校正后 p 值 < 0.01
        abs(avg_log2FC) > 0.5      # FC绝对值 > 0.5（约1.4倍变化）
      ) %>%  
      arrange(p_val_adj) %>%       # 按 p_val_adj 升序排序（p值越小越靠前）
      head(n_genes)                # 取前 n_genes 个最显著的基因
    
    # 然后筛选在两组中表达比例都≥50%的基因
    DotPlot(
      subset_data,
      features = cell_deg$gene,
      group.by = "group",
      cols = c("blue", "red"),
      dot.scale = 6
    ) +
      theme(
        axis.text.x = element_text(angle = 45, hjust = 1, face = "italic"),
        legend.position = "right"
      ) +
      labs(x = "Genes", y = "Group") +
      scale_y_discrete(limits = c("Healthy", "ASS-ILD")) +
      ggtitle(paste("Top", n_genes, "Most different DEGs (p_adj < 0.01)"))
  }
  
  # 修改报告生成函数中的点图调用(唯一需要修改的地方)
  generate_deg_report <- function(cell_type) {
    cell_deg <- combined_deg %>% filter(cell_type == !!cell_type)
    
    if(nrow(cell_deg) > 0) {
      # 创建可视化组合(将plot_dotplot改为plot_dotplot_top_genes)
      p1 <- plot_volcano(cell_deg, paste(cell_type, "ASS-ILD vs Healthy"))
      p2 <- plot_dotplot_top_genes(combined_seurat, cell_type)  # 改为新的点图函数
      
      # 调整布局
      combined_plot <- p1 / p2 + 
        plot_layout(heights = c(1, 1.2))
      
      # 保存结果
      ggsave(
        paste0("DEG_results/", cell_type, "_DEG_report.pdf"),
        combined_plot,
        width = 12,
        height = 14
      )
    }
  }
  # 其余代码保持不变...
  
  # 生成报告时强制遍历所有目标细胞类型
  for(cell in target_celltypes) {  # 修改为遍历指定细胞类型
    if(cell %in% unique(combined_deg$cell_type)) {  # 只在有结果时生成
      generate_deg_report(cell)
    } else {
      message(paste("No DEGs found for", cell))
    }
  }
}
#组间差异分析对两个组所有的基因
{
  # 确认分组信息已正确添加
  table(combined_seurat$group)  # 应显示ASS-ILD和Healthy组的细胞数量
  Idents(combined_seurat) <- "group"  # 设置分组为当前细胞标识
  # 使用Wilcoxon秩和检验（推荐用于单细胞数据）
  deg_markers <- FindMarkers(
    object = combined_seurat,
    ident.1 = "ASS-ILD",      # 实验组
    ident.2 = "Healthy",      # 对照组
    test.use = "wilcox",      # 非参数检验适合稀疏数据
    logfc.threshold = 0.25,   # 最小log2倍变化阈值[6](@ref)
    min.pct = 0.25,            # 基因至少在10%的细胞中表达[4](@ref)
    only.pos = FALSE          # 同时输出上调/下调基因
  )
  
  # 添加基因名称列
  deg_markers$gene <- rownames(deg_markers)
  
  # 按调整p值和fold change排序
  deg_markers <- deg_markers[order(deg_markers$p_val_adj, -abs(deg_markers$avg_log2FC)), ]
  
  # 保存结果
  write.csv(deg_markers, "Group_DEG_Results.csv", row.names = FALSE)
  # 定义显著性阈值（可调整）
  
  sig_genes3 <- subset(deg_markers, 
                       p_val_adj < 0.01& 
                         abs(avg_log2FC) > 1)
  
                
                
  
}
#GO/KEEG(富集分析)
{
  {
    #富集分析
    # 加载必要的包
    library(clusterProfiler)
    library(org.Hs.eg.db)
    library(enrichplot)
    library(ggplot2)
    library(dplyr)
    library(Seurat)
    
    # 创建富集分析主目录
    if(!dir.exists("Enrichment_Results")) dir.create("Enrichment_Results")
    
    # 修改后的富集分析函数（移除热图）
    run_enrichment_analysis <- function(deg_df, analysis_type, comparison, 
                                        p_val_cutoff = 0.05, logfc_cutoff = 0.25) {
      # 参数校验
      if(!all(c("gene", "avg_log2FC", "p_val_adj") %in% colnames(deg_df))) {
        stop("Input dataframe must contain columns: gene, avg_log2FC, p_val_adj")
      }
      
      # 创建专属目录
      output_dir <- file.path("Enrichment_Results", analysis_type)
      if(!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)
      
      # 基因筛选
      sig_genes <- deg_df %>%
        filter(p_val_adj < p_val_cutoff & abs(avg_log2FC) > logfc_cutoff) %>%
        pull(gene)
      
      if(length(sig_genes) < 5) {
        message("Insufficient genes for enrichment in ", analysis_type)
        return()
      }
      
      # ID转换
      gene.df <- tryCatch({
        bitr(sig_genes, 
             fromType = "SYMBOL",
             toType = "ENTREZID",
             OrgDb = org.Hs.eg.db)
      }, error = function(e) {
        message("ID conversion failed for ", analysis_type, ": ", e$message)
        return(NULL)
      })
      
      if(is.null(gene.df) || nrow(gene.df) < 5) {
        message("Insufficient mapped genes for ", analysis_type)
        return()
      }
      
      # GO富集分析
      run_go_enrichment <- function(ontology) {
        enrichGO(gene = gene.df$SYMBOL,
                 OrgDb = org.Hs.eg.db,
                 keyType = "SYMBOL",
                 ont = ontology,
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.05,
                 qvalueCutoff = 0.2,
                 readable = TRUE)
      }
      
      go_results <- list(
        BP = run_go_enrichment("BP"),
        MF = run_go_enrichment("MF"),
        CC = run_go_enrichment("CC")
      )
      
      # KEGG通路分析
      kegg_result <- tryCatch({
        enrichKEGG(
          gene = gene.df$ENTREZID,
          organism = "hsa",
          pvalueCutoff = 0.05,
          pAdjustMethod = "BH",
          qvalueCutoff = 0.2
        )
      }, error = function(e) {
        message("KEGG analysis failed for ", analysis_type, ": ", e$message)
        return(NULL)
      })
      
      # 保存结果
      save_enrichment_results <- function(result, prefix) {
        if(!is.null(result) && nrow(result) > 0) {
          write.csv(as.data.frame(result), 
                    file.path(output_dir, paste0(prefix, "_enrichment.csv")))
        }
      }
      
      save_enrichment_results(go_results$BP, "GO_BP")
      save_enrichment_results(go_results$MF, "GO_MF")
      save_enrichment_results(go_results$CC, "GO_CC")
      save_enrichment_results(kegg_result, "KEGG")
      
      # 修改后的可视化函数（移除热图）
      generate_enrichment_plots <- function(result, title_suffix, type) {
        if(!is.null(result) && nrow(result) > 0) {
          try({
            # 生成显示名称
            display_name <- ifelse(
              grepl("CellType_", analysis_type),
              gsub("CellType_", "", analysis_type),
              "Group Comparison"
            )
            
            # 点图
            p1 <- dotplot(result, showCategory=15) + 
              ggtitle(paste(title_suffix, "in", display_name)) +
              theme(plot.title = element_text(hjust = 0.5))
            
            # 网络图
            p2 <- emapplot(pairwise_termsim(result), showCategory = 15) +
              ggtitle(paste(title_suffix, "Network in", display_name))
            
            # 保存图形
            ggsave(file.path(output_dir, paste0(type, "_dotplot.pdf")), 
                   p1, width=10, height=8)
            ggsave(file.path(output_dir, paste0(type, "_network.pdf")), 
                   p2, width=12, height=10)
          })
        }
      }
      
      # 生成可视化
      generate_enrichment_plots(go_results$BP, "Biological Process", "GO_BP")
      generate_enrichment_plots(go_results$MF, "Molecular Function", "GO_MF")
      generate_enrichment_plots(go_results$CC, "Cellular Component", "GO_CC")
      generate_enrichment_plots(kegg_result, "KEGG Pathways", "KEGG")
      
      # KEGG通路可视化
      if(!is.null(kegg_result) && nrow(kegg_result) > 0) {
        try({
          gene_fc <- setNames(deg_df$avg_log2FC[match(gene.df$SYMBOL, deg_df$gene)],
                              gene.df$ENTREZID)
          
          for(i in 1:min(3, nrow(kegg_result))) {
            pathview(
              gene.data = gene_fc,
              pathway.id = kegg_result$ID[i],
              species = "hsa",
              gene.idtype = "ENTREZID",
              limit = list(gene=2, cpd=1),
              out.suffix = paste0(analysis_type, "_pathway"),
              kegg.dir = output_dir
            )
            
            # 移动生成的图形文件
            png_pattern <- paste0("hsa", kegg_result$ID[i], ".*", analysis_type, "_pathway.png")
            png_files <- list.files(pattern = png_pattern)
            file.rename(png_files, 
                        file.path(output_dir, paste0(kegg_result$ID[i], "_pathway.png")))
          }
        })
      }
    }
    
    # 主分析流程
    # 1. 细胞类型特异性分析
    if(exists("combined_deg")) {
      all_cell_types <- unique(combined_deg$cell_type)
      for(cell in all_cell_types) {
        current_deg <- combined_deg %>% filter(cell_type == cell)
        if(nrow(current_deg) > 0) {
          run_enrichment_analysis(
            current_deg,
            analysis_type = paste0("CellType_", cell),
            comparison = "ASS-ILD_vs_Healthy",
            p_val_cutoff = 0.05,
            logfc_cutoff = 0.25
          )
        }
      }
    }
    
    # 2. 组间整体分析
    if(exists("deg_markers") && nrow(deg_markers) > 0) {
      run_enrichment_analysis(
        deg_markers,
        analysis_type = "Group_Comparison",
        comparison = "ASS-ILD_vs_Healthy",
        p_val_cutoff = 0.01,
        logfc_cutoff = 0.25
      )
    }
  }
}
#细胞通讯(总体)
{
  if(!exists("combined_seurat") || !"SingleR_labels" %in% colnames(combined_seurat@meta.data)) {
    stop("请先完成细胞注释部分（SingleR注释）")
  }
  # 准备输入数据
  {
    # 提取细胞类型标签（使用SingleR注释结果）
    cell_labels <- combined_seurat$SingleR_labels
    
    # 提取标准化后的表达矩阵（使用log归一化数据）
    data.input <- GetAssayData(combined_seurat, assay = "RNA", layer = "data")  # Seurat v5使用layer参数
    
    # 创建CellChat对象（使用稀疏矩阵节省内存）
    cellchat <- createCellChat(
      object = data.input,
      meta = data.frame(
        labels = cell_labels,
        group = combined_seurat$group,  # 保留分组信息
        row.names = colnames(data.input)
      ),
      group.by = "labels",
      do.sparse = TRUE
    )
    
    # 清理临时变量释放内存
    rm(data.input); gc()
  }
  # 设置CellChat数据库
  {
    # 使用人类数据库（如果是小鼠数据则用CellChatDB.mouse）
    CellChatDB <- CellChatDB.human  
    
    # 可选：仅保留分泌信号（根据研究需求调整）
    # CellChatDB <- subsetDB(CellChatDB, search = "Secreted Signaling")
    
    # 将数据库链接到CellChat对象
    cellchat@DB <- CellChatDB
  }
  # 预处理和通讯网络推断
  {
    # 子集数据（仅保留在CellChatDB中定义的配体-受体）
    cellchat <- subsetData(cellchat) 
    
    # 识别过表达基因和相互作用
    cellchat <- identifyOverExpressedGenes(cellchat)
    cellchat <- identifyOverExpressedInteractions(cellchat)
    
    # 计算通讯概率（使用trimmed mean减少异常值影响）
    cellchat <- computeCommunProb(cellchat, type = "truncatedMean", trim = 0.1)
    
    # 过滤低可信度通讯（至少10个细胞）
    cellchat <- filterCommunication(cellchat, min.cells = 10)
    
    # 计算通路水平的通讯网络
    cellchat <- computeCommunProbPathway(cellchat)
    
    # 整合通讯网络
    cellchat <- aggregateNet(cellchat)
  }
  groupSize <- as.numeric(table(cellchat@idents))
  par(mfrow = c(1,2))
  netVisual_circle(cellchat@net$count, vertex.weight = groupSize, title.name = "Interaction Counts")
  netVisual_circle(cellchat@net$weight, vertex.weight = groupSize, title.name = "Interaction Strength")
  # 2. 信号通路可视化（示例使用MIF通路）
  pathways.show <- "MIF"
  vertex.receiver <- c(1,3)  # 指定接收细胞群的索引
  
  # 层次图
  netVisual_aggregate(cellchat, 
                      signaling = pathways.show, 
                      vertex.receiver = vertex.receiver,
                      layout = "hierarchy")
  
  # 弦图
  netVisual_aggregate(cellchat, 
                      signaling = pathways.show,
                      layout = "chord")
  
  # 热图
  netVisual_heatmap(cellchat, 
                    signaling = pathways.show, 
                    color.heatmap = "Reds")
  
  # 3. 配体-受体贡献度分析
  netAnalysis_contribution(cellchat, signaling = pathways.show)
  
  # 4. 基因表达可视化
  plotGeneExpression(cellchat, signaling = pathways.show)
  ### 系统分析
  # 1. 网络中心性分析
  cellchat <- netAnalysis_computeCentrality(cellchat)
  netAnalysis_signalingRole_network(cellchat, signaling = pathways.show)
  
  {
    # 加载包
    library(NMF)
    library(ggalluvial)
    library(doParallel)
    options(repr.plot.height = 8, repr.plot.width = 8)
    selectK(cellchat, pattern = "incoming")
    nPatterns = 4
    cellchat <- identifyCommunicationPatterns(cellchat, 
                                              pattern = "incoming", 
                                              k = nPatterns,
                                              width = 8, height = 18)
    netAnalysis_river(cellchat, pattern = "incoming")
    netAnalysis_dot(cellchat, pattern = "incoming")
    
  }
}
save.image("细胞通讯整体.RData")
saveRDS(cellchat,file = "细胞通讯整体.rds")
cellchat <-  readRDS("细胞通讯整体.rds")
#细胞通讯分别创建cellchat对象
{
  # 细胞通讯差异分析
  {
    # 检查是否已完成细胞注释和分组
    if(!exists("combined_seurat") || 
       !"SingleR_labels" %in% colnames(combined_seurat@meta.data) ||
       !"group" %in% colnames(combined_seurat@meta.data)) {
      stop("请先完成细胞注释和分组分析部分")
    }
    
    # 提取分组信息（根据您的代码1，分组为ASS-ILD和Healthy）
    meta_healthy <- subset(combined_seurat@meta.data, group == "Healthy")
    meta_disease <- subset(combined_seurat@meta.data, group == "ASS-ILD")  # 修改为您的疾病组名称
    
    # 提取标准化表达数据
    data.input <- GetAssayData(combined_seurat, assay = "RNA", layer = "data")  # Seurat v5使用layer参数
    
    # 创建分组的数据子集
    data_healthy <- data.input[, rownames(meta_healthy)]
    data_disease <- data.input[, rownames(meta_disease)]
    
    # 创建分组的CellChat对象
    cellchat_healthy <- createCellChat(
      data_healthy, 
      meta = meta_healthy, 
      group.by = "SingleR_labels",
      do.sparse = TRUE  # 使用稀疏矩阵节省内存
    )
    
    cellchat_disease <- createCellChat(
      data_disease, 
      meta = meta_disease, 
      group.by = "SingleR_labels",
      do.sparse = TRUE
    )
    
    # 统一使用人类数据库
    CellChatDB <- CellChatDB.human  # 使用完整数据库
    cellchat_healthy@DB <- CellChatDB
    cellchat_disease@DB <- CellChatDB
    
    # 定义预处理函数（添加错误处理和内存管理）
    process_group <- function(cellchat_obj) {
      tryCatch({
        cellchat_obj <- subsetData(cellchat_obj)
        cellchat_obj <- identifyOverExpressedGenes(cellchat_obj)
        cellchat_obj <- identifyOverExpressedInteractions(cellchat_obj)
        
        # 使用trimmed mean减少异常值影响
        cellchat_obj <- computeCommunProb(cellchat_obj, type = "truncatedMean", trim = 0.1)
        
        # 过滤低可信度通讯
        cellchat_obj <- filterCommunication(cellchat_obj, min.cells = 10)
        
        # 计算通路水平通讯
        cellchat_obj <- computeCommunProbPathway(cellchat_obj)
        
        # 整合网络
        cellchat_obj <- aggregateNet(cellchat_obj)
        
        return(cellchat_obj)
      }, error = function(e) {
        message("处理过程中出错: ", e$message)
        return(NULL)
      })
    }
    
    # 处理各组数据
    cellchat_healthy <- process_group(cellchat_healthy)
    cellchat_disease <- process_group(cellchat_disease)
  }
  saveRDS(cellchat_healthy,file = "cellchat_healthy.rds")
  saveRDS(cellchat_disease,file = "cellchat_disease.rds")
}
cellchat_healthy <-  readRDS("cellchat_healthy.rds")
cellchat_disease <- readRDS("cellchat_disease.rds")
#组间强度比较
{
  # 创建分组对象列表
  cellchat.list <- list(Healthy = cellchat_healthy, Disease = cellchat_disease)
  # 合并对象（关键步骤）
  cellchat_merged <- mergeCellChat(
    cellchat.list, 
    add.names = names(cellchat.list),
    cell.prefix = TRUE # 保留原始细胞ID前缀
  )
  # 通讯数量对比（柱状图）
  p1 <- compareInteractions(cellchat_merged, 
                            measure = "count",
                            color.use = c("#1f78b4", "red"),
                            group = c("Healthy", "Disease"))
  
  # 通讯强度对比（柱状图）
  p2 <- compareInteractions(cellchat_merged,
                            measure = "weight",
                            color.use = c("#1f78b4", "red"),
                            group = c("Healthy", "Disease"))
  
  # 组合图形
  p1 | p2
}
#细胞类型
{
  # 检查健康组细胞类型
  cat("健康组细胞类型:\n")
  print(table(cellchat_healthy@idents))
  
  # 检查疾病组细胞类型
  cat("\n疾病组细胞类型:\n")
  print(table(cellchat_disease@idents))
  
  # 可视化细胞类型组成差异
  celltype_comparison <- merge(
    as.data.frame(table(cellchat_healthy@idents), responseName = "Healthy"),
    as.data.frame(table(cellchat_disease@idents), responseName = "Disease"),
    by = "Var1",  # 按细胞类型合并
    all = TRUE    # 保留所有细胞类型
  )
  colnames(celltype_comparison)[1] <- "CellType"
  
  # 替换NA为0
  celltype_comparison[is.na(celltype_comparison)] <- 0
  
  # 绘制堆叠条形图
  library(reshape2)
  library(ggplot2)
  melted_data <- melt(celltype_comparison, id.vars = "CellType")
  ggplot(melted_data, aes(x = variable, y = value, fill = CellType)) +
    geom_bar(stat = "identity", position = "fill") +
    labs(x = "Group", y = "Proportion", title = "健康组与疾病组细胞类型组成比较") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
  
  # 找出健康组特有的细胞类型
  healthy_unique <- setdiff(levels(cellchat_healthy@idents), levels(cellchat_disease@idents))
  cat("\n健康组特有细胞类型:", if(length(healthy_unique)>0) healthy_unique else "无", "\n")
  
  # 找出疾病组特有的细胞类型
  disease_unique <- setdiff(levels(cellchat_disease@idents), levels(cellchat_healthy@idents))
  cat("疾病组特有细胞类型:", if(length(disease_unique)>0) disease_unique else "无", "\n")
  
  # 找出共有细胞类型
  common_types <- intersect(levels(cellchat_healthy@idents), levels(cellchat_disease@idents))
  cat("\n共有细胞类型:", common_types, "\n")
  
  # 计算比例变化
  celltype_comparison$ratio <- celltype_comparison$Disease / celltype_comparison$Healthy
  celltype_comparison$ratio[is.infinite(celltype_comparison$ratio)] <- NA
  
  # 找出在疾病组中比例显著增加的细胞类型（top 3）
  increased <- celltype_comparison[order(-celltype_comparison$ratio, na.last = TRUE), ]
  cat("\n疾病组中比例增加最多的细胞类型（Top 3）:\n")
  print(head(increased[, c("CellType", "Healthy", "Disease", "ratio")], 3))
  
  # 找出在疾病组中比例显著减少的细胞类型（top 3）
  decreased <- celltype_comparison[order(celltype_comparison$ratio, na.last = TRUE), ]
  cat("\n疾病组中比例减少最多的细胞类型（Top 3）:\n")
  print(head(decreased[, c("CellType", "Healthy", "Disease", "ratio")], 3))
  
  # 可选：绘制比例变化点图
  ggplot(celltype_comparison, aes(x = Healthy, y = Disease, color = CellType)) +
    geom_point(size = 3) +
    geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
    scale_x_log10() + scale_y_log10() +
    labs(title = "细胞类型数量变化", x = "健康组细胞数量", y = "疾病组细胞数量") +
    theme_minimal()
}
#特有细胞通讯
{
  # 健康组特有细胞分析
  if(length(healthy_unique) > 0) {
    cat("\n分析健康组特有细胞类型的通讯模式...\n")
    
    # 提取健康组特有细胞的通讯网络
    net_healthy_source <- subsetCommunication(cellchat_healthy, 
                                              sources.use = healthy_unique,
                                              targets.use = levels(cellchat_healthy@idents))
    
    net_healthy_target <- subsetCommunication(cellchat_healthy, 
                                              sources.use = levels(cellchat_healthy@idents),
                                              targets.use = healthy_unique)
    
    # 可视化特有细胞通讯
    for(ct in healthy_unique) {
      # 作为源的通讯可视化
      if(nrow(subset(net_healthy_source, source == ct)) > 0) {
        netVisual_chord_gene(cellchat_healthy, 
                             sources.use = ct, 
                             targets.use = setdiff(levels(cellchat_healthy@idents), ct),
                             lab.cex = 0.8, 
                             legend.pos.y = 30,
                             title.name = paste("健康组特有细胞:", ct, "作为源的通讯"))
      } else {
        cat("注意：健康组特有细胞", ct, "作为源无显著通讯\n")
      }
      
      # 作为靶的通讯可视化
      if(nrow(subset(net_healthy_target, target == ct)) > 0) {
        netVisual_chord_gene(cellchat_healthy, 
                             sources.use = setdiff(levels(cellchat_healthy@idents), ct), 
                             targets.use = ct,
                             lab.cex = 0.8, 
                             legend.pos.y = 30,
                             title.name = paste("健康组特有细胞:", ct, "作为靶的通讯"))
      } else {
        cat("注意：健康组特有细胞", ct, "作为靶无显著通讯\n")
      }
    }
    
    # 保存结果
    write.csv(net_healthy_source, "healthy_unique_as_source.csv")
    write.csv(net_healthy_target, "healthy_unique_as_target.csv")
  } else {
    cat("\n健康组无特有细胞类型\n")
  }
  
  # 疾病组特有细胞分析
  if(length(disease_unique) > 0) {
    cat("\n分析疾病组特有细胞类型的通讯模式...\n")
    
    # 提取疾病组特有细胞的通讯网络
    net_disease_source <- subsetCommunication(cellchat_disease, 
                                              sources.use = disease_unique,
                                              targets.use = levels(cellchat_disease@idents))
    
    net_disease_target <- subsetCommunication(cellchat_disease, 
                                              sources.use = levels(cellchat_disease@idents),
                                              targets.use = disease_unique)
    
    # 可视化特有细胞通讯
    for(ct in disease_unique) {
      # 作为源的通讯可视化
      if(nrow(subset(net_disease_source, source == ct)) > 0) {
        netVisual_chord_gene(cellchat_disease, 
                             sources.use = ct, 
                             targets.use = setdiff(levels(cellchat_disease@idents), ct),
                             lab.cex = 0.8, 
                             legend.pos.y = 30,
                             title.name = paste("疾病组特有细胞:", ct, "作为源的通讯"))
      } else {
        cat("注意：疾病组特有细胞", ct, "作为源无显著通讯\n")
      }
      
      # 作为靶的通讯可视化
      if(nrow(subset(net_disease_target, target == ct)) > 0) {
        netVisual_chord_gene(cellchat_disease, 
                             sources.use = setdiff(levels(cellchat_disease@idents), ct), 
                             targets.use = ct,
                             lab.cex = 0.8, 
                             legend.pos.y = 30,
                             title.name = paste("疾病组特有细胞:", ct, "作为靶的通讯"))
      } else {
        cat("注意：疾病组特有细胞", ct, "作为靶无显著通讯\n")
      }
    }
    
    # 保存结果
    write.csv(net_disease_source, "disease_unique_as_source.csv")
    write.csv(net_disease_target, "disease_unique_as_target.csv")
  } else {
    cat("\n疾病组无特有细胞类型\n")
  }
}
#通路比较
{
  # 1. 获取各组显著通路
  pathways_healthy <- cellchat_healthy@netP$pathways
  pathways_disease <- cellchat_disease@netP$pathways
  
  # 打印各组通路信息
  cat("\n健康组通路 (", length(pathways_healthy), "条):\n", paste(pathways_healthy, collapse=", "))
  cat("\n疾病组通路 (", length(pathways_disease), "条):\n", paste(pathways_disease, collapse=", "))
  
  # 2. 系统性通路比较分析
  pathway_analysis <- list(
    # 特有通路
    unique_healthy = setdiff(pathways_healthy, pathways_disease),
    unique_disease = setdiff(pathways_disease, pathways_healthy),
    
    # 共有通路
    common = intersect(pathways_healthy, pathways_disease)
  )
  
  # 3. 格式化输出结果
  cat("\n\n=== 通路分析结果 ===\n")
  for (category in names(pathway_analysis)) {
    cat("\n", gsub("_", " ", toupper(category)), " (", 
        length(pathway_analysis[[category]]), "条):\n")
    if (length(pathway_analysis[[category]]) > 0) {
      cat(paste(pathway_analysis[[category]], collapse=", "), "\n")
    } else {
      cat("无\n")
    }
  }
  
  # 4. 可视化特有通路
  # 健康组特有通路可视化
  if (length(pathway_analysis$unique_healthy) > 0) {
    pdf("healthy_unique_pathways.pdf", width = 10, height = 8)
    for (pathway in pathway_analysis$unique_healthy) {
      tryCatch({
        write(pathway)
        netVisual_aggregate(cellchat_healthy, 
                            signaling = pathway, 
                            layout = "circle",
        )
        title(main = paste("healthy only:", pathway), 
              line = -1, outer = TRUE)
      }, error = function(e) {
        cat("可视化通路", pathway, "时出错:", e$message, "\n")
      })
    }
    dev.off()
    cat("\n健康组特有通路可视化已保存到: healthy_unique_pathways.pdf\n")
  }
  
  # 疾病组特有通路可视化
  if (length(pathway_analysis$unique_disease) > 0) {
    pdf("disease_unique_pathways.pdf", width = 10, height = 8)
    for (pathway in pathway_analysis$unique_disease) {
      tryCatch({
        netVisual_aggregate(cellchat_disease, 
                            signaling = pathway, 
                            layout = "circle",
                            title = paste("疾病组特有通路:", pathway))
        title(main = paste("disease only:", pathway), 
              line = -1, outer = TRUE)
      }, error = function(e) {
        cat("可视化通路", pathway, "时出错:", e$message, "\n")
      })
    }
    dev.off()
    cat("疾病组特有通路可视化已保存到: disease_unique_pathways.pdf\n")
  }else{
    cat("无通路")
  }
  
  # 5. 共有通路活性比较
  if (length(pathway_analysis$common) > 0) {
    # 创建通路活性数据框
    pathway_activity <- data.frame(
      Pathway = pathway_analysis$common,
      Healthy = sapply(pathway_analysis$common, function(p) {
        sum(cellchat_healthy@netP$prob[,,p], na.rm = TRUE)
      }),
      Disease = sapply(pathway_analysis$common, function(p) {
        sum(cellchat_disease@netP$prob[,,p], na.rm = TRUE)
      }),
      stringsAsFactors = FALSE
    )
    
    # 计算差异和变化倍数
    pathway_activity$Difference <- pathway_activity$Disease - pathway_activity$Healthy
    pathway_activity$FoldChange <- ifelse(pathway_activity$Healthy > 0,
                                          pathway_activity$Disease / pathway_activity$Healthy,
                                          NA)
    
    # 筛选显著差异通路（基于标准差）
    diff_threshold <- mean(abs(pathway_activity$Difference), na.rm = TRUE) + 
      sd(abs(pathway_activity$Difference), na.rm = TRUE)
    sig_pathways <- pathway_activity[abs(pathway_activity$Difference) > diff_threshold, ]
    
    # 可视化显著通路
    if (nrow(sig_pathways) > 0) {
      # 热图可视化
      library(ComplexHeatmap)
      mat <- as.matrix(sig_pathways[, c("Healthy", "Disease")])
      rownames(mat) <- sig_pathways$Pathway
      
      pdf("common_pathways_heatmap.pdf", width = 8, height = 10)
      draw(Heatmap(mat, 
                   name = "Activity",
                   col = colorRamp2(c(min(mat), median(mat), max(mat)), 
                                    c("blue", "white", "red")),
                   column_title = "共有通路活性比较",
                   row_names_gp = gpar(fontsize = 8),
                   column_names_rot = 45))
      dev.off()
      
      # 差异条形图
      library(ggplot2)
      sig_pathways$Direction <- ifelse(sig_pathways$Difference > 0, "Up", "Down")
      ggplot(sig_pathways, aes(x = reorder(Pathway, Difference), y = Difference, fill = Direction)) +
        geom_bar(stat = "identity") +
        coord_flip() +
        labs(x = "Pathway", y = "Disease - Healthy", 
             title = "共有通路活性差异") +
        theme_minimal() +
        scale_fill_manual(values = c("Up" = "red", "Down" = "blue"))
      ggsave("common_pathways_difference.pdf", width = 8, height = 6)
    }
    
    # 保存结果
    write.csv(pathway_activity, "pathway_activity_comparison.csv", row.names = FALSE)
    cat("\n共有通路活性比较结果已保存到: pathway_activity_comparison.csv\n")
  } else {
    cat("\n警告：未发现共有通路，跳过活性比较分析\n")
  }
  
  # 6. 通路网络差异可视化（示例展示共有通路）
  if (length(pathway_analysis$common) > 0) {
    common_pathways <- head(pathway_analysis$common, 46)
    pdf("pathway_network_comparison.pdf", width = 14, height = 6)
    par(mfrow = c(1, 2))
    for (pathway in common_pathways) {
      tryCatch({
        netVisual_aggregate(cellchat_healthy, 
                            signaling = pathway, 
                            layout = "circle"
        )
        title(main = paste("healthy-VS-disease:", pathway), line = -1, outer = TRUE)
        netVisual_aggregate(cellchat_disease, 
                            signaling = pathway, 
                            layout = "circle")
        
      }, error = function(e) {
        cat("可视化通路", pathway, "时出错:", e$message, "\n")
      })
    }
    dev.off()
    cat("\n通路网络对比图已保存到: pathway_network_comparison.pdf\n")
  }
}
#配受体比较
{
  # 创建结果目录
  results_dir <- "CellChat_Healthy_vs_Disease_Results"
  if (!dir.exists(results_dir)) {
    dir.create(results_dir)
  }
  
  # 提取健康组和疾病组的配体-受体通信数据
  lr_healthy <- subsetCommunication(cellchat_healthy)
  lr_disease <- subsetCommunication(cellchat_disease)
  
  # 检查数据有效性
  if(!is.null(lr_healthy) && !is.null(lr_disease) && 
     nrow(lr_healthy) > 0 && nrow(lr_disease) > 0) {
    
    # 添加组别标识
    lr_healthy$group <- "Healthy"
    lr_disease$group <- "Disease"
    
    # 合并数据
    lr_combined <- rbind(lr_healthy, lr_disease)
    
    # 计算组间差异和变异指标
    lr_diff <- lr_combined %>%
      group_by(source, target, pathway_name, ligand, receptor) %>%
      summarise(
        Healthy_prob = mean(prob[group == "Healthy"], na.rm = TRUE),
        Disease_prob = mean(prob[group == "Disease"], na.rm = TRUE),
        # 组间差异
        Difference = Disease_prob - Healthy_prob,
        # 变化倍数（防止除零错误）
        FoldChange = ifelse(Healthy_prob > 0, 
                            Disease_prob / Healthy_prob,
                            NA),
        # 变异指标
        Variability = abs(Difference),  # 绝对差异值
        .groups = 'drop'
      ) %>%
      arrange(desc(Variability))  # 按变异程度排序
    
    # 筛选显著差异的配体-受体对（前10%）
    sig_threshold <- quantile(lr_diff$Variability, 0.9, na.rm = TRUE)
    sig_lr_diff <- lr_diff %>% 
      filter(Variability >= sig_threshold)
    
    cat("\n显著差异的配体-受体对数量:", nrow(sig_lr_diff), "\n")
    
    # 可视化分析（仅当存在显著差异时）
    if(nrow(sig_lr_diff) > 0) {
      # 取变异最大的前20对（或全部如果不足20）
      top_lr <- head(sig_lr_diff, min(20, nrow(sig_lr_diff)))
      # 1. 两组通讯概率对比条形图
      library(dplyr)
      library(tidyr)  # pivot_longer需要这个包
      
      # 修改后的代码
      plot_data_prob <- top_lr %>%
        dplyr::select(ligand, receptor, Healthy_prob, Disease_prob) %>%
        pivot_longer(cols = c(Healthy_prob, Disease_prob),
                     names_to = "Group", 
                     values_to = "Probability") %>%
        mutate(Group = factor(Group, 
                              levels = c("Healthy_prob", "Disease_prob"),
                              labels = c("Healthy", "Disease")),
               LR_pair = paste0(ligand, "-", receptor))
      
      p_prob <- ggplot(plot_data_prob, 
                       aes(x = reorder(LR_pair, Probability), 
                           y = Probability, 
                           fill = Group)) +
        geom_bar(stat = "identity", position = position_dodge(), width = 0.7) +
        coord_flip() +
        scale_fill_manual(values = c("Healthy" = "#1b9e77", 
                                     "Disease" = "#d95f02")) +
        labs(x = "Ligand-Receptor Pair", 
             y = "Communication Probability",
             title = "Top Differential LR Pairs: Communication Probability") +
        theme_minimal() +
        theme(legend.position = "top",
              axis.text.y = element_text(size = 8))
      
      # 2. 差异程度点图（按通路着色）
      p_diff <- ggplot(top_lr, 
                       aes(x = reorder(paste0(ligand, "-", receptor), Difference), 
                           y = Difference,
                           color = pathway_name,
                           size = abs(Difference))) +
        geom_point() +
        geom_hline(yintercept = 0, linetype = "dashed") +
        coord_flip() +
        scale_color_brewer(palette = "Set2", name = "Pathway") +
        scale_size_continuous(name = "Absolute Difference") +
        labs(x = "Ligand-Receptor Pair", 
             y = "Probability Difference (Disease - Healthy)",
             title = "Top Differential LR Pairs: Communication Differences") +
        theme_minimal() +
        theme(axis.text.y = element_text(size = 8),
              legend.position = "right")
      
      # 3. 变化倍数热图
      plot_fc <- top_lr %>%
        mutate(LR_pair = paste0(ligand, "-", receptor),
               log2FC = log2(FoldChange))  # 对数转换
      
      p_fc <- ggplot(plot_fc, 
                     aes(x = "Comparison", 
                         y = reorder(LR_pair, log2FC), 
                         fill = log2FC)) +
        geom_tile() +
        scale_fill_gradient2(low = "blue", mid = "white", high = "red",
                             midpoint = 0, name = "log2(Fold Change)") +
        labs(x = "", y = "Ligand-Receptor Pair",
             title = "Top Differential LR Pairs: Fold Change") +
        theme_minimal() +
        theme(axis.text.x = element_blank(),
              axis.text.y = element_text(size = 8))
      
      # 保存所有图形
      ggsave(file.path(results_dir, "1_Probability_Comparison.pdf"), p_prob, 
             width = 10, height = 8)
      ggsave(file.path(results_dir, "2_Difference_Dotplot.pdf"), p_diff, 
             width = 12, height = 8)
      ggsave(file.path(results_dir, "3_FoldChange_Heatmap.pdf"), p_fc, 
             width = 8, height = 8)
      
      
      # 保存结果数据
      write.csv(lr_diff, file.path(results_dir, "All_LR_Pair_Differences.csv"), 
                row.names = FALSE)
      write.csv(sig_lr_diff, file.path(results_dir, "Significant_LR_Pair_Differences.csv"), 
                row.names = FALSE)
      
      # 输出结果说明
      cat("\n分析结果已保存到:", results_dir, "\n")
      cat("包含以下文件:\n")
      cat("- 1_Probability_Comparison.pdf: 两组通讯概率对比条形图\n")
      cat("- 2_Difference_Dotplot.pdf: 差异程度点图（按通路着色）\n")
      cat("- 3_FoldChange_Heatmap.pdf: 变化倍数热图\n")
      cat("- All_LR_Pair_Differences.csv: 所有配体-受体对差异数据\n")
      cat("- Significant_LR_Pair_Differences.csv: 显著差异配体-受体对数据\n")
      
      # 额外保存每组特有的LR对
      healthy_only <- anti_join(lr_healthy, lr_disease, 
                                by = c("ligand", "receptor", "pathway_name"))
      disease_only <- anti_join(lr_disease, lr_healthy, 
                                by = c("ligand", "receptor", "pathway_name"))
      
      if(nrow(healthy_only) > 0) {
        write.csv(healthy_only, file.path(results_dir, "Healthy_Specific_LR_Pairs.csv"), 
                  row.names = FALSE)
        cat("- Healthy_Specific_LR_Pairs.csv: 健康组特有的配体-受体对\n")
      }
      if(nrow(disease_only) > 0) {
        write.csv(disease_only, file.path(results_dir, "Disease_Specific_LR_Pairs.csv"), 
                  row.names = FALSE)
        cat("- Disease_Specific_LR_Pairs.csv: 疾病组特有的配体-受体对\n")
      }
    }
  } else {
    cat("\n警告：配体-受体对数据不完整，无法进行比较分析。\n")
    if(is.null(lr_healthy) || nrow(lr_healthy) == 0) {
      cat("健康组数据存在问题\n")
    }
    if(is.null(lr_disease) || nrow(lr_disease) == 0) {
      cat("疾病组数据存在问题\n")
    }
  }
}
#信号流交流分析
{
  signaling_dir <- file.path(results_dir, "Signaling_Flow_Analysis")
  # 确保结果目录存在
  if (!dir.exists(results_dir)) {
    dir.create(results_dir, recursive = TRUE)
  }
  if (!dir.exists(signaling_dir)) {
    dir.create(signaling_dir, recursive = TRUE)
  }
  
  # 计算网络中心性指标
  cellchat_healthy <- netAnalysis_computeCentrality(cellchat_healthy)
  cellchat_disease <- netAnalysis_computeCentrality(cellchat_disease)
  
  # 获取共有通路
  common_pathways <- intersect(cellchat_healthy@netP$pathways, 
                               cellchat_disease@netP$pathways)
  
  if (length(common_pathways) > 0) {
    cat("\n开始分析", length(common_pathways), "个共有信号通路...\n")
    
    # 为每个通路创建子目录
    for (p in common_pathways) {
      # 清理通路名称用于文件名
      clean_pathway_name <- gsub("[^[:alnum:]_]", "_", p)
      pathway_dir <- file.path(signaling_dir, clean_pathway_name)
      
      # 确保通路目录存在
      if (!dir.exists(pathway_dir)) {
        dir.create(pathway_dir, recursive = TRUE)
      }
      
      # 生成PDF文件名
      pdf_file <- file.path(pathway_dir, paste0("Signaling_Flow_", clean_pathway_name, ".pdf"))
      
      # 创建PDF文件
      pdf_success <- tryCatch({
        while (dev.cur() > 1) { dev.off() } # 关闭已打开的图形设备
        pdf(pdf_file, width = 12, height = 8)
        TRUE
      }, error = function(e) {
        message("无法创建PDF文件 ", pdf_file, ": ", e$message)
        FALSE
      })
      
      if (pdf_success) {
        tryCatch({
          # 健康组信号流
          p1 <- try({
            plot_obj <- netAnalysis_signalingRole_network(cellchat_healthy, 
                                                          signaling = p,
                                                          width = 20,
                                                          height = 12)
            if (!inherits(plot_obj, "try-error")) {
              grid.text(paste("Healthy -", p), x = 0.5, y = 0.95, 
                        gp = gpar(fontsize = 16, fontface = "bold"))
              print(plot_obj)
            }
          })
          
          # 疾病组信号流
          p2 <- try({
            plot_obj <- netAnalysis_signalingRole_network(cellchat_disease, 
                                                          signaling = p,
                                                          width = 20,
                                                          height = 12)
            if (!inherits(plot_obj, "try-error")) {
              grid.text(paste("Disease -", p), x = 0.5, y = 0.95, 
                        gp = gpar(fontsize = 16, fontface = "bold"))
              print(plot_obj)
            }
          })
          
          # 添加分页
          grid.newpage()
        }, finally = {
          if (dev.cur() > 1) dev.off()
        })
      }
      
      # 保存PNG格式图片
      try({
        png(file.path(pathway_dir, paste0("Healthy_", clean_pathway_name, ".png")), 
            width = 1200, height = 800)
        plot_obj <- netAnalysis_signalingRole_network(cellchat_healthy, signaling = p)
        grid.text(paste("Healthy -", p), x = 0.5, y = 0.95, 
                  gp = gpar(fontsize = 16, fontface = "bold"))
        print(plot_obj)
        dev.off()
      }, silent = TRUE)
      
      try({
        png(file.path(pathway_dir, paste0("Disease_", clean_pathway_name, ".png")), 
            width = 1200, height = 800)
        plot_obj <- netAnalysis_signalingRole_network(cellchat_disease, signaling = p)
        grid.text(paste("Disease -", p), x = 0.5, y = 0.95, 
                  gp = gpar(fontsize = 16, fontface = "bold"))
        print(plot_obj)
        dev.off()
      }, silent = TRUE)
    }
    
    # 生成分析报告
    try({
      sink(file.path(signaling_dir, "Analysis_Report.txt"))
      cat("信号流比较分析报告\n")
      cat("生成日期:", date(), "\n\n")
      cat("分析组别: Healthy vs Disease\n")
      cat("分析的共有通路数量:", length(common_pathways), "\n")
      cat("通路列表:\n")
      print(common_pathways)
      cat("\n可视化结果包含:\n")
      cat("- PDF文件包含两组的信号流对比图\n")
      cat("- PNG文件包含各组的独立信号流图\n")
      sink()
    }, silent = TRUE)
    
    cat("\n分析完成! 结果保存在:", signaling_dir, "\n")
  } else {
    cat("\n警告: 没有共有通路，跳过信号流分析。\n")
  }
}
#细胞通讯组间差异
{
  # 加载所需包
  library(CellChat)
  library(patchwork)
  library(ComplexHeatmap)
  
  # 1. 数据准备和合并 =======================================================
  
  # 找出所有细胞类型的并集
  all_celltypes <- unique(c(
    levels(cellchat_healthy@idents),
    levels(cellchat_disease@idents)
  ))
  
  # 提升各组的细胞标签到统一标准
  cellchat_healthy <- liftCellChat(cellchat_healthy, group.new = all_celltypes)
  cellchat_disease <- liftCellChat(cellchat_disease, group.new = all_celltypes)
  
  # 合并对象
  object.list <- list(
    Healthy = cellchat_healthy,
    Disease = cellchat_disease
  )
  
  cellchat.merged <- mergeCellChat(
    object.list,
    add.names = names(object.list),
    cell.prefix = TRUE
  )
  
  # 2. 比较细胞通讯网络的一般特征 ==========================================
  
  # 比较交互总数和交互强度（合并展示）
  gg1 <- compareInteractions(cellchat.merged, show.legend = F, group = 1:2, size.text = 14)
  gg2 <- compareInteractions(cellchat.merged, show.legend = F, group = 1:2, measure = "weight", size.text = 14)
  combined_plot <- gg1 + gg2 + plot_annotation(title = "Interaction Counts and Weights Comparison")
  print(combined_plot)
  
  # 网络图比较（并排展示）
  weight.max <- getMaxWeight(object.list, attribute = c("idents", "count"))
  par(mfrow = c(1,2), xpd = TRUE, mar = c(4,4,4,4))
  for (i in 1:length(object.list)) {
    netVisual_circle(object.list[[i]]@net$count, 
                     vertex.weight = as.numeric(table(object.list[[i]]@idents)),
                     weight.scale = T, 
                     edge.weight.max = weight.max[2], 
                     edge.width.max = 10, 
                     title.name = paste0(names(object.list)[i], "\nInteraction Counts"))
  }
  
  # 差异网络可视化
  par(mfrow = c(1,1), xpd=TRUE)
  netVisual_diffInteraction(cellchat.merged, comparison = c(1,2), 
                            weight.scale = T, title.name = "Disease vs Healthy")
  
  # 差异网络热图
  heatmap_plot <- netVisual_heatmap(cellchat.merged, comparison = c(1,2), 
                                    title.name = "Disease vs Healthy", 
                                    font.size = 12)
  print(heatmap_plot)
  
  # 3. 信号网络中心性分析 =================================================
  
  # 计算网络中心性
  for (i in 1:length(object.list)) {
    object.list[[i]] <- netAnalysis_computeCentrality(object.list[[i]])
  }
  
  # 信号角色散点图比较
  weight.MinMax <- c(
    min(sapply(object.list, function(x) sum(x@net$count))),
    max(sapply(object.list, function(x) sum(x@net$count)))
  )
  
  scatter_plots <- list()
  for (i in 1:length(object.list)) {
    scatter_plots[[i]] <- netAnalysis_signalingRole_scatter(
      object.list[[i]],
      title = names(object.list)[i],
      weight.MinMax = weight.MinMax
    )
  }
  wrap_plots(scatter_plots, ncol = 2)
  
  # 4. 通路活性比较 =======================================================
  
  # 比较各通路的整体信息流
  rank_plot <- rankNet(cellchat.merged, mode = "comparison", 
                       comparison = c(1,2), stacked = T,
                       font.size = 12, title = "Pathway Activity Comparison")
  print(rank_plot)
  
  # 5. 共有通路可视化 =====================================================
  # 获取共有通路
  common_pathways <- intersect(
    cellchat_healthy@netP$pathways,
    cellchat_disease@netP$pathways
  )
  
  cat("Number of common pathways:", length(common_pathways), "\n")
  
  # 批量保存共有通路气泡图
  if(length(common_pathways) > 0) {
    for(pathway in common_pathways){
      # 创建组合图形
      combined_bubble <- netVisual_bubble(cellchat_healthy, signaling = pathway, 
                                          title = paste("Healthy -", pathway)) + 
        netVisual_bubble(cellchat_disease, signaling = pathway, 
                         title = paste("Disease -", pathway))
      
      # 保存PDF
      pdf(file = paste0("CommonPathway_", pathway, ".pdf"), width = 12, height = 6)
      print(combined_bubble)
      dev.off()
    }
  }
  
  # 6. 重点通路详细分析 ===================================================
  if(length(common_pathways) > 0) {
    pathway.show <- common_pathways[1]
    
    # 比较两组中的通路活性 - 圆形图
    par(mfrow = c(1,2), xpd=TRUE, mar = c(0, 0, 2, 0))
    for (i in 1:length(object.list)) {
      netVisual_aggregate(object.list[[i]], signaling = pathway.show, 
                          layout = "circle", edge.weight.max = weight.max[1], 
                          edge.width.max = 10, 
                          title = paste(names(object.list)[i], "\n", pathway.show))
    }
    
    # 和弦图可视化比较
    par(mfrow = c(1,2), xpd=TRUE, mar = c(0, 0, 2, 0))
    for (i in 1:length(object.list)) {
      netVisual_aggregate(object.list[[i]], signaling = pathway.show, 
                          layout = "chord", 
                          title.name = paste(names(object.list)[i], "\n", pathway.show),
                          small.gap = 1,
                          big.gap = 5,
                          reduce = -1)
    }
    
    # 比较信号基因表达
    gene_plot <- plotGeneExpression(cellchat.merged, signaling = pathway.show, 
                                    split.by = "datasets", colors.ggplot = TRUE)
    gene_plot <- gene_plot + ggtitle(paste("Gene Expression -", pathway.show))
    print(gene_plot)
  }
  
  # 7. 结果汇总输出 =======================================================
  cat("\nAnalysis completed!\n")
  cat("Key comparisons between Healthy and Disease groups:\n")
  cat("- Total interactions:", sum(cellchat_healthy@net$count), "(Healthy) vs", 
      sum(cellchat_disease@net$count), "(Disease)\n")
  cat("- Number of significant pathways:", length(common_pathways), "\n")
  if(length(common_pathways) > 0) {
    cat("- Example pathway analysis performed for:", pathway.show, "\n")
  }
}